###---###---###---###---###---###---###---###---###---###---
# Kerice Doten-Snitker
# 
# Dissertation chapter: Diffusion of Expulsions
# Snippet: working on path distances along routes
# 
# Built under R 3.6.1 
# Platform: x86_64-apple-darwin15.6.0 (64-bit)
###---###---###---###---###---###---###---###---###---###---

# set locale to preempt UTF-8 issues with the mydata
#Sys.getlocale()
#Sys.setlocale("LC_CTYPE", "en_US.UTF-8")

#for windows
#Sys.setlocale(category = "LC_ALL", locale = "English_United States.1252")

# add packages with the pacman library, which installs if not installed
library(pacman)
# use the here package for improving replicability
p_load(here)
# for data manipulation
p_load(magrittr, plyr, tidyverse, reshape2)
# for geographical analyses and calculations
p_load(geosphere, raster, gdistance, rgdal, rgeos)
p_load(igraph,riverdist)

set.seed(1418)


# resource on routing: 
# Rowlingson, Barry. "Line to Graph to Route." http://rpubs.com/geospacedman/routing

### Import data with the cities and their coordinates -----------------------

# import csv of settlement coordinates (easier for some calculations)
points_coords_meters <- read.csv(here::here("Diffusion", "out_files", "points_coords_meters.csv"), header=TRUE, row.names=1, stringsAsFactors=FALSE)

# shp of settlement coordinates (easier for some calculations)
cities <- readOGR(here::here("GIS_data","Vectors", "GderJuden_cleaned"), "GdJ_cleaned_nocov_meters", GDAL1_integer64_policy = TRUE, stringsAsFactors=FALSE)

###---###---###---###---###---###---###---###---###---###---
# Distances between cities using medieval pilgrimage routes ----
###---###---###---###---###---###---###---###---###---###---

# The routes used here were shared by Bossak and Welford as the individual pilgrimage routes. I cleaned them by snapping segment ends together and merging them all into a single lines shapefile.

### Workflow in QGIS 3.4 Madeira ###---###---###---###---###---###

# QGIS has dialog box options to snap points to lines, and to calculate the distances between the old points and the new snapped points, but it's not easily replicable.
# 1) snap cities to road and river networks with Processing Toolbox ->  Vector Geometry Tools -> "Snap geometries to layer" : tolerance of 100 km, "prefer closest point, do not add vertices"
# 2) Either save in the snapping dialog, or create temp layer and then save
# 3) Processing -> Network Analysis -> Point to Layer (one at a time)

# Because snapping distances have to be calculated individually with the QGIS dialog (well, they could be done all at once using Python within QGIS, but...), it's just much easier to use an R-based process that will do all these steps.

#### Using R: with riverdist package ------------------------------
crs_EPSG3035 <- "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
# set the tolerance 1 because the unit is meters
roads_net <- line2network(here::here("GIS_data","Vectors", "medieval_routes"), layer="routes_all-segments_addvert_meters", tolerance=1, reproject=crs_EPSG3035)
# do some built-in processing - y/n questions
roads_net_cleaned <- cleanup(roads_net)
# dissolve:y split:y insert verts:y min dist:1000 
# mouth segment:118 mouth vertex: 78 accept:y remove:n remove:n braiding:n
saveRDS(roads_net_cleaned, here::here("in_data", "roads_net_cleaned_meters.rds"))
roads_net_lim <- trimtopoints(x=points_coords_meters$XCOORD, y=points_coords_meters$YCOORD, rivers=roads_net_cleaned, method="buffer", dist=10000)

# "snap" the cities to locations on the road network
# points_road <- xy2segvert(x=points_coords_meters$XCOORD, y=points_coords_meters$YCOORD, rivers=roads_net_cleaned)
# this one will keep settlement metadata attached, like name
# setts_road <- pointshp2segvert(path=here::here("GIS_data","Vectors", "GderJuden_cleaned"), layer="GdJ_cleaned_nocov_meters",
#                                      rivers=roads_net_lim)
# again with shp, but already snapped with QGIS
setts_snap <- pointshp2segvert(path=here::here("GIS_data","Vectors", "GderJuden_cleaned"), 
                               layer="GdJ_cleaned_snapped_road100km", rivers=roads_net_lim)


# check what it looks like
zoomtoseg(seg=c(32, 43), rivers=roads_net_lim)
points(points_coords$XCOORD, points_coords$YCOORD, pch=16, col="red")
riverpoints(seg=setts_snap$seg, vert=setts_snap$vert, rivers=roads_net_lim, pch=15, 
            col="blue")

dist_rnet_matrix <- riverdistancemat(setts_snap$seg, setts_snap$vert,
  roads_net_lim, ID=setts_snap$PlaceName, 
  stopiferror=FALSE, algorithm="Dijkstra")

write.csv(dist_rnet_matrix,here::here("Diffusion", "out_files", "distances_meters_roads.csv"), row.names=TRUE)

#Utrecht to Muenster
setts_snap$seg[setts_snap$PlaceName=="Utrecht"]
setts_snap$vert[setts_snap$PlaceName=="Utrecht"]
setts_snap$seg[setts_snap$PlaceName=="Muenster"]
setts_snap$vert[setts_snap$PlaceName=="Muenster"]
riverdistance(startseg=23, startvert=12, endseg=32, endvert=16, rivers=roads_net_lim, map=TRUE)

# calc distances btwn cities and locations snapped to roads

setts_on <- readOGR(here::here("GIS_data","Vectors", "GderJuden_cleaned"), "GdJ_cleaned_snapped_road100km", GDAL1_integer64_policy = TRUE)
setts_on_df <- data.frame(setts_on)
# duplicate Oberehnheim rows...need to fix this upstream
# as.integer(rownames(setts_on_df[setts_on_df$PlaceName=="Oberehnheim",])) #317
# setts_on_df <- setts_on_df[-318,] # choose 7.498, 48.454

dist2snap_roads <- dist(setts_on_df[,c(3,4)], setts_on_df[,c(5,6)])/1000 # distances in meters, so divide for km
dist2snap_roads <- as.data.frame(dist2snap_roads)
rownames(dist2snap_roads) <- setts_on_df$PlaceName

#write.csv(dist2snap_roads,here::here("Diffusion", "out_files", "dist2snap_meters_roads.csv"), row.names=TRUE)

###---###---###---###---###---###---###---###---###---###---###---###
### Route-based distances between cities: add rivers ----------------
###---###---###---###---###---###---###---###---###---###---###---###

# The majority of cities are nearby BW routes, but some are visibly along alternatives: rivers! (many routes were following rivers). I augment the route data with river shapefiles from Natural Earth II.

### in QGIS ###---###---###---###---###---###

# QGIS has dialog box options to snap points to lines, and to calculate the distances between the old points and the new snapped points, but it's not easily replicable.
# 1) snap cities to road and river networks with Processing Toolbox ->  Vector Geometry Tools -> "Snap geometries to layer" : tolerance of 100 km, "prefer closest point, do not add vertices"
# 2) Either save in the snapping dialog, or create temp layer and then save
setts_onboth <- readOGR(here::here("GIS_data","Vectors", "GderJuden_cleaned"), "GdJ_cleaned_snapped_mainriv01road100km", GDAL1_integer64_policy = TRUE)
# 3) Processing -> Network Analysis -> Point to Layer (one at a time)
# 4) multiply resulting cost by 110 to transform to kms
# for later: get distance between original and updated lat-long

#### R: with riverdist ------------------------------

# this will take a few minutes to load - 8 min
# perhaps should have higher tolerance, since the underlying file was simplified
tstart <- Sys.time()
rivroad_net <- line2network(here::here("GIS_data","Vectors", "routes_plus_main-rivers"), layer="routes_plus_main-rivers-01_addvert_meters", tolerance=10, reproject=crs_EPSG3035)
tend <- Sys.time()
saveRDS(rivroad_net, here::here("Diffusion","out_files", "mainriv01road_meters_net.rds"))
# do some built-in processing - y/n questions
tstart <- Sys.time() # takes >2hrs
rivroad_net_cleaned <- cleanup(rivroad_net)
# dissolve:y split:y insert verts:y min dist:1000 
# mouth seg:391 mouth vertex:78 accept:y remove:n remove:n braiding:n
#saveRDS(rivroad_net_cleaned, here::here("in_data", "rivroad_meters_net_cleaned.rds"))
# 26 regions need manual review - add connections for most
rivroad_net_lim <- trimtopoints(x=points_coords_meters$XCOORD, y=points_coords_meters$YCOORD, rivers=rivroad_net_cleaned, method="buffer", dist=20000) #already clipped

# saving the network file seems to strip its projection info, as when I work from the .rds to snap points below it tells me that the projections are different and the points will be reprojected before snapping...
saveRDS(rivroad_net_lim, here::here("Diffusion","out_files", "mainriv01road_meters_net_lim.rds"))

# "snap" the settlements to locations on the road network
#points_rivroad <- xy2segvert(x=points_coords_meters$XCOORD, y=points_coords_meters$YCOORD, rivers=rivroad_net_lim)
# this one will keep settlement metadata attached, like name
setts_rivroad <- pointshp2segvert(path=here::here("GIS_data","Vectors", "GderJuden_cleaned"), layer="GdJ_cleaned_nocov_meters",
                               rivers=rivroad_net_lim)
# again with shp, but already snapped with QGIS
setts_snap_rr <- pointshp2segvert(path=here::here("GIS_data","Vectors", "GderJuden_cleaned"), layer="GdJ_cleaned_snapped_mainriv01road100km", 
                                  rivers=rivroad_net_lim)

# check what it looks like
zoomtoseg(seg=c(32, 43), rivers=rivroad_net_lim)
points(points_coords_meters$XCOORD, points_coords_meters$YCOORD, pch=16, col="red")
riverpoints(seg=setts_snap_rr$seg, vert=setts_snap_rr$vert, rivers=rivroad_net_lim, pch=15, 
            col="blue")

tstart <- Sys.time()
dist_rivroad_matrix <- riverdistancemat(setts_snap_rr$seg, setts_snap_rr$vert,
                                        rivroad_net_lim, ID=setts_snap_rr$PlaceName, 
                                        stopiferror=FALSE, algorithm="Dijkstra")
tend <- Sys.time()

# check a route
setts_snap_rr$seg[setts_snap_rr$PlaceName=="Duisburg"]
setts_snap_rr$vert[setts_snap_rr$PlaceName=="Duisburg"]
setts_snap_rr$seg[setts_snap_rr$PlaceName=="Heilbronn"]
setts_snap_rr$vert[setts_snap_rr$PlaceName=="Heilbronn"]
riverdistance(startseg=250, startvert=182, endseg=208, endvert=99, rivers=rivroad_net_lim, map=TRUE)

write.csv(dist_rivroad_matrix,here::here("Diffusion", "out_files", "distances_mainriv01road_meters.csv"), row.names=TRUE)

# calc distances btwn cities and locations snapped to roads
setts_onboth$XCOORD <- cities$coords.x1
setts_onboth$YCOORD <- cities$coords.x2
setts_onboth_df <- data.frame(setts_onboth)

dist2snap_rivroad <- pointDistance(setts_onboth_df[,c(3,4)], # the original coords
                          setts_onboth_df[,c(5,6)], 
                          lonlat=FALSE, allpairs=FALSE)/1000 # the snapped coords; distances in meters, so divide for km
dist2snap_rivroad <- as.data.frame(dist2snap_rivroad)
rownames(dist2snap_rivroad) <- setts_onboth_df$PlaceName

#write.csv(dist2snap_rivroad,here::here("Diffusion", "out_files", "dist2snap_mainriv01road_rivroad.csv"), row.names=TRUE)
